#
# BAM_LAT_PARTIES.r
#
library(rjags)
library(foreign)
library(coda)
#
LAT <- read.dta("/Users/ssaiegh/Mac/Representation Gap/Paper PA/stimuli.dta",convert.factors=FALSE)
attach(LAT)
#
LA_Self <- self
H.Chavez <- chavez
A.Uribe <- uribe
C.Kirchner <- cristina
D.Scioli <- scioli
E.Morales <- evo
Lula <- lula
D.Rousseff <- dilma
J.Serra <- serra
S.Pinera <- pinera
M.Bachelet <- bachelet
J.M.Santos <- jm_santos
A.Mockus <- mockus
F.Calderon <- calderon
A.Garcia <- alan
O.Humala <- ollanta
J.Mujica <- mujica
L.A.Lacalle <- lacalle
B.Obama <- obama
R.Zapatero <- zapatero
# 
FPV <- fpv
PJ <- pj
UCR <- ucr
PRO <- pro
PT <- pt
PMDB <- pmdb
PSDB <- psdb
DEM <- dem
FA <- fa
PNAC <- pn
PCOL <-pc
PDC <-pdc
PPD <- ppd
PS <- ps
RN <- rn
UDI <- udi
PRI <- pri
PAN <- pan
PRD <- prd
MAS <- mas
PPB <- ppb
MIR <- mir
PdeU <- p_de_u
PCC <- pcc
PLC <- plc
APRA <- pap
G.PERU <- gana_peru
#
LA_TT <- cbind(LA_Self,H.Chavez,A.Uribe,C.Kirchner,D.Scioli,E.Morales,Lula,D.Rousseff,
J.Serra,S.Pinera,M.Bachelet,J.M.Santos,A.Mockus,F.Calderon,A.Garcia,
O.Humala,J.Mujica,L.A.Lacalle,B.Obama,R.Zapatero,
FPV,PJ,UCR,PRO,PT,PMDB,PSDB,DEM,FA,PNAC,PCOL,PDC,PPD,PS,RN,UDI,PRI,PAN,PRD,
MAS,PPB,MIR,PdeU,PCC,PLC,APRA,G.PERU)
#
# Replace missing entries with "NA":
LA_TT[LA_TT==98 | LA_TT==99] <- NA
LA_TT <- sweep(LA_TT,1,rowMeans(LA_TT,na.rm=T),"-")
# Set the minimum number of placements that a respondent needs to
# provide to be included in the scaling (default = 3):
LA_cutoff <- 3
LA_cutoff <- ncol(LA_TT) - LA_cutoff
#
LA_TT <- LA_TT[rowSums(is.na(LA_TT))<=LA_cutoff,]
#
self <- LA_TT[,1]
LA_TT <- LA_TT[,-1]
#
N <- nrow(LA_TT)
q <- ncol(LA_TT)
z <- LA_TT
#
alpha.starts <- rnorm(N,0,0.3)
beta.starts <- rnorm(N,1,0.3)
#
zhat.starts <- rnorm(q,0,0.3)
zhat.starts[1] <- runif(1,-1.1,-0.9)
zhat.starts[2] <- runif(1,0.9,1.1)
#
inits <- function() {list (zhat=zhat.starts, a=alpha.starts, b=beta.starts)}
#
BAM.sim <- jags.model('/Users/ssaiegh/Mac/Representation Gap/Paper PA/BAM_JAGScode.bug',
	data = list('z' = z,'q' = q,'N' = N),inits = inits, n.chains = 2, n.adapt = 10000) 
update(BAM.sim, n.iter = 10000)
#
zhat <- coda.samples(BAM.sim, 'zhat', 2500, thin=5)
a <- coda.samples(BAM.sim, 'a', 2500, thin=5)
b <- coda.samples(BAM.sim, 'b', 2500, thin=5)
tauj <- coda.samples(BAM.sim, 'tauj', 2500, thin=5)
#
#
#  STIMULI ESTIMATES
#
summary(zhat)
#
#
#  DIAGNOSTICS
#
plot(zhat[,1:3])
geweke.diag(zhat)
gelman.diag(zhat)
#
#
